The superfluid fountain effect in a Bose-Einstein condensate 
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We consider a simple experimental setup, based on a harmonic confinement, where a Bose-Einstein 
condensate and a thermal cloud of weakly interacting alkali atoms are trapped in two different vessels 
connected by a narrow channel. Using the classical field approximation, as described in J. Phys. 
B 40, Rl (2007) and optimized in Phys. Rev. A 81, 013629 (2010) for an arbitrary trapping 
potential, we theoretically investigate the analog of the celebrated superfluid helium fountain effect. 
We show that this thermo-mechanical effect might indeed be observed in this system. By analyzing 
the dynamics of the system, we are able to identify the superfluid and normal components of the 
flow as well as to distinguish the condensate fraction from the superfluid component. We show that 
the superfluid component can easily flow from the colder vessel to the hotter one while the normal 
component is practically blocked in the latter. 



I. INTRODUCTION 

The experimental discovery of superfluidity in helium 
II by Kapitsa [T] and Allen and Misener [5] in 1938 has 
triggered a great theoretical interest in this phenomenon. 
One of the most spectacular effects related to superflu- 
idity of helium II is its ability to flow through narrow 
channels with apparently zero viscosity. Extensive stud- 
ies of this system were very important for the founda- 
tion of the theory of Bose and Fermi quantum liquids. 
In this system however, even at the lowest temperatures, 
the strong interactions between the helium atoms deplete 
the population of the Bose-Einstein condensate to about 
10% of the total mass whereas the superfluid fraction is 
almost 100%. 

The situation is substantially different with dilute ul- 
tracold atomic gases. The first implementation of a 
Bose-Einstein condensation [31 Sj in alkali atoms has 
opened new possibilities to explore Bose quantum liq- 
uids at much higher level of control. Indeed, contrary 
to liquid helium, large condensate fractions are routinely 
obtained with dilute atomic gases as the atoms are very 
weakly interacting. To date, many phenomena previ- 
ously observed in liquid helium below the lambda point 
have found their experimental counterpart with ultracold 
alkali gases even if the experimental evidence of superflu- 
idity in atomic condensates has been a very challenging 
task. One of the main signatures of superfluid flow is the 
generation of quantized vortices when the system is set 
into rotation. After many efforts such quantized vortices, 
and also arrays of vortices, were observed in atomic con- 
densates |5rt^ . Observation of the first sound [3 |9] , of 
scissor modes [10] or of the critical velocity [11] beyond 
which the superfluid flow breaks down, are other exam- 



ples of the manifestation of this spectacular macroscopic 
quantum phenomenon in trapped ultracold atomic sys- 
tems. 

In addition to the above-mentioned properties, helium 
II exhibits also a very unusual feature related to the flow 
of heat. Variations of temperature propagate in this sys- 
tem in a form of waves known as the second sound. Both 
these extraordinary features, i.e. non viscous flow and 
unusual heat transport, manifest themselves in full glory 
in the helium fountain effect, called also the thermo- 
mechanical effect. Its first observation was reported by 
Allen and Jones [12]. In their original setup, the lower 
part of a U-tube packed with fine emery powder was im- 
mersed into a vessel containing liquid helium II. A tem- 
perature gradient was created by shining a light beam 
on the powder which got heated due to light absorption. 
As a result of the temperature gradient, a superfluid flow 
is generated from the cold liquid helium reservoir to the 
hotter region. This flow can be so strong that a jet of he- 
lium is forced up through the vertical part of the U-tube 
to a height of several centimeters, hence the fountain ef- 
fect name. 

Up to now, there exists many different experimental 
implementations of this spectacular effect and one of 
them is shown in FigjlJ A small vessel, connected to a 
bulb filled with emery powder forming a very fine capil- 
lary net, is immersed in the container of liquid helium II. 
When the electric heater is off the superfluid liquid flows 
freely through the capillary net in the bulb and fills in the 
small vessel as shown in panel (A) of Figjl] If now the su- 
perfluid helium inside the small vessel is heated then the 
level of the liquid in the vessel increases above the level 
of the liquid in the big container, see panel (B). A con- 
tinuous heating sustains the flow from the colder part of 
the system to the hotter one, an observation at variance 
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with our ordinary everyday life experience. Eventually 
liquid helium reaches the top of the small vessel where it 
forms the helium fountain, see panel (C). 

The explanation for this counter-intuitive thermo- 
mechanical effect is closely related to the notion of the 
second sound and to the two-fluid model developed by 
Tisza and Landau [131 [H] • This approach assumes the 
existence of two co-existing components of the liquid he- 
lium: the superfluid and the normal one. The normal 
component is viscous and can transport heat. On con- 
trary, the superfluid component has no viscosity and can- 
not transport heat. Because it is viscous, the normal 
component cannot flow through the capillary net but the 
superfluid can. Heat transport is thus forbidden because 
it can only be carried by the normal component. As a 
consequence, the system cannot reach thermal equilib- 
rium and the temperature in the reservoir keeps smaller 
than the temperature in the small vessel. Only mechani- 
cal equilibrium is allowed, i.e. the chemical potentials in 
both vessels have to equilibrate. Heating of the superfluid 
component inside the small vessel leads to a reduction of 
the chemical potential per unit mass in this vessel and 
this reduction is compensated by the flow of the super- 
fluid component from the reservoir. 

The two-fluid model for helium H assumes a local ther- 
mal equilibrium which signifies a hydrodynamic regime 
where the collision time is the shortest time scale. If 
this is indeed the case for superfluid helium H, which 
is a strongly interacting system, it is generally not for 
trapped ultracold dilute atomic gases where reaching this 
regime proves extremely difficult. For example, second 
sound has only been observed recently [TBI. As a con- 
sequence, the usual two-fluid model fails to apply and, 
up to our knowledge, there are no theoretical predictions 
about heat transport in weakly-interacting atomic con- 
densates. It is not even obvious if an effect similar to the 
helium fountain can be observed in these dilute atomic 
systems. This is probably why the analog of this classic 
textbook phenomenon has not yet been studied in the 
context of ultracold trapped atomic gases. 

Nevertheless the question of the nature of heat trans- 
port in these weakly-interacting atomic condensates 
seems to be well posed. There are not many experi- 
ments where a non-equilibrium transfer of atoms related 
to temperature differences have been studied. We should 
recall here the experiment of the MIT group, where dis- 
tillation of a condensate was observed [T^. The authors 
studied how the superfluid system 'discovers' the exis- 
tence of a dynamically-created global minimum of the 
trapping potential and how the system gets to this mini- 
mum. Theoretical studies of the corresponding ID situa- 
tion suggested different dynamical behaviors of the ther- 
mal fraction and of the superfluid component which, in 
some sense, resemble the fountain effect |18) . 

In the present work we theoretically study the non- 
equilibrium dynamics of a Bose-Einstein condensate 
which is driven by a temperature gradient. We will show 
that an effect qualitatively very similar to the helium 



fountain may be observed in experiments with trapped 
ultracold dilute atomic gases. 




FIG. 1. A cartoon picture showing the idea of the superfluid 
fountain experiment. 



The paper is organized as follows: In Sec. [TTj we briefly 
introduce the classical field approximation. Sec. [Hi] de- 
scribes the system under consideration and our numerical 
procedure. Then in Sec. |IV| we present and analyze our 
numerical data. We show that the thermo-mechanical ef- 
fect is indeed present in our system and we highlight the 
importance of distinguishing between the superfluid, nor- 
mal, condensate and thermal components of the system. 
Finally, we give in Sec. IVlsomc concluding remarks. 



II. CLASSICAL FIELDS APPROXIMATION 

There exist different effective methods to describe and 
study dynamical effects in condensates at nonzero tem- 
perature. For example, the Zaremba-Nikuni-Griffin for- 
malism assumes a splitting of the system into a conden- 
sate and a thermal cloud [19j whereas different versions of 
the classical fields method describe both the condensate 
and the thermal cloud by a single Gross-Pitaevskii equa- 
tion [5DH23]. Here, we will use the classical fields method 
as described in |25j and optimized in |26j for arbitrary 
trapping potentials. Since this latter paper describes 
quite extensively how to prepare the classical fields in a 
thermodynamic equilibrium state characterized by tem- 
perature T, total number of atoms N, condensate frac- 
tion Nf)/N {Nq being the number of condensed atoms), 
and scattering length a, we will just briefly describe here 
the main ingredients of this approach. 



A. Formalism 

We start with the usual bosonic field operator \E'(r,t) 
which annihilates an atom at point r and time t. It obeys 
the standard bosonic commutation relations: 



^(r,i),^+(r',i) =(5(r-r') 

*+(r,i),^+(r',t)]=0 

[*(r,t),*(r',i)]=0, (1) 
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and evolves according to the Heisenberg equation of mo- 
tion: 



2m 

+g*+(r,i)*(r,<)*(r,i) 



*(r,i) 



(2) 



where Vtr{r,t) is the (possibly) time-dependent trapping 
potential and g = AT:h?a/m is the coupling constant ex- 
pressed in terms of the s-wave scattering length a. 

The field operator itself can be expanded in a basis of 
one-particle wave functions 0Q(r), where a denotes the 
set of all necessary one-particle quantum numbers: 



(3) 



In the presence of a trap, a natural choice for the 
one-particle modes (jia would be the harmonic oscillator 
modes, otherwise one generally uses plane wave states. 
The classical fields method is an extension of the Bogoli- 
ubov idea to finite temperatures and gives some micro- 
scopic basis to the two-fluid model. The main idea is to 
assume that modes (j)a in expansion ([s]) having an energy 
Ea less than a certain cut-off energy Ec are macroscop- 
ically occupied and, consequently, to replace all corre- 
sponding annihilation operators by c-number amplitudes: 

*(r,t)~ ^ 0„(r)a„(i)+ ^ 0„(r)a„(i) . (4) 



-Eq>-Ec 



Assuming further that the second sum in Q is small and 
can be neglected, the field operator 'I'(r,t) is turned into 
a classical complex wave function: 

i>{v,t)^^{Y,t)= J2 Mr)aa.it). (5) 

In this way, both the condensate and a thermal cloud of 
atoms, interacting with each other, will be described by 
a single classical field "^{r,t). Injecting ([5| into ([2]), we 
obtain the equation of motion for the classical field: 



Zm 



+g**(r,t)*(r,i)*(r,t) 



(6) 



In numerical implementations, one controls a total en- 
ergy, a number of macroscopically occupied modes 0q 
and a value of gN. The energy-truncation constraint 
Ea < Ec is usually implemented by solving Eq. ([6| on a 
rectangular grid using the Fast Fourier Transform tech- 
nique. The spatial grid step determines the maximal mo- 
mentum per particle, and hence the energy, in the system 
whereas the use of the Fourier transform implies projec- 
tion in momentum space. 

Equation ^ looks identical to the usual Gross- 
Pitaevskii equation describing a Bose-Einstein conden- 
sate at zero temperature. However, the interpretation of 



the complex wave function \E'(r, t) is here different. It de- 
scribes all the atoms in the system. Therefore, the ques- 
tion arises on how to extract all these modes out of the 
time-evolving classical field ^(r, t). For this purpose, we 
follow the definition of Bose-Einstein condensation orig- 
inally proposed by Penrose and Onsager [28] where the 
condensate is assigned to be described by the eigenvec- 
tor corresponding to the dominant eigenvalue of the one- 
particle density matrix. This one-particle density matrix 
reads: 



eW(r,r';t) = ^vl/(r,i)xp*(r',i), 



1 



(7) 



and obviously corresponds to a pure state with all atoms 
in the condensate mode. This is because Eq. ^ is the 
spectral decomposition of the one-particle density ma- 
trix. To extract the modes out of the classical field some 
kind of averaging is needed. In a typical experiment, one 
generally measures the column density integrated along 
some direction. We implement here the same type of 
procedure and define the coarse-grained density matrix: 



eix,y,x',y';t) 



1 

N 



dz^{x,y,z,t)^*{x',y',z,t), 



(8) 

from which we extract the corresponding eigenvalues in 
order to apply the Penrose-Onsager criterion. We have 
tested the ability of this averaging procedure for a clas- 
sical field at thermal equilibrium in a harmonic trap by 
comparing it to the results obtained when the original 
density matrix is averaged long enough over time. With 
a 1% accuracy, both methods give the same results. 

Solving the eigenvalue problem for the coarse-grained 
density matrix ([s]) leads to the decomposition: 

K 

Q{x,y,x,y;t) = ^nk{t) (pk{x,y,t) (pl{x' ,y' ,t) , (9) 

k=Q 

where the relative occupation numbers nk{t) — Nklt)/N 
of the orthonormal macroscopically occupied modes 
are ordered according to no(t) > ni{t) >(...)> nxit). 
For future convenience, we define the eigenmodes of 
the coarse-grained one-particle density matrix which are 
normalized to the relative occupation numbers of these 
modes and the corresponding one-particle density matrix 
Pk- 



^k{x,y,t) = 



^k{x,y,t) 



Qk{x,y,x',y';t) = tpk{x,y,t)ipl{x' ,y' ,t) , (10) 



such that g — X^s^o ^^'^ = X]fe=i 8k, the conden- 
sate being described by Qq. 

According to the standard definition, the condensate 
wave function corresponds to fc = and the thermal den- 
sity is simply: 



PT{x,y,t) = g{x,y,x,y;t) ~ |-0o(a;, y, t)|^ 



(11) 



4 



In an equilibrium situation, the relative occupation num- 
bers rifc do not depend on time. In this case, the to- 
tal number of atoms is determined from the smallest 
eigenvalue of the one-particle density matrix through 
ukN = ricut, where Ucut ~ 0.46 for the 3D harmonic 
oscillator [29], from which one can infer the value of the 
interaction strength g. The temperature T of the sys- 
tem is then given by the energy of this highest occupied 
mode. 

Let us note that any initial state evolving with the 
Gross-Pitaevskii equation reaches a state of thermal equi- 
librium characterized by a temperature, a total number 
of particles and an interaction strength g. However to 
obtain an equilibrium classical field for a given set of pa- 
rameters, one has to properly choose the energy of the 
initial state and the cut-off parameter Ec- This task is 
time consuming because the temperature and the num- 
ber of particles can only be assigned to the field after 
the equilibrium is reached. To speed up the prepara- 
tion of the initial equilibrium state, we first solve the 
self-consistent Hartree-Fock model (SCHFM) [35]. This 
allows us to estimate quite accurately the energy of the 
state for a given temperature and particle number. The 
detailed description of this procedure can be found in 
[26] . Here we only recall the SCHFM equations: 

Po{r)^-[ti-Vtr{v)-2gpT{v)] (12) 
9 

fir, p) - (e[pV2m+V,(r)-p]/fe«T _ (^3) 

PT(r) = ^ .93/2 (14) 

Veir)^Vtr{r) + 2gpo{r) + 2gpTir) (15) 
p^gpo{0) + 2g pt{0) + Vtr{0) , (16) 

where 

is the thermal de Broglie wavelength. The 53/2 (z) func- 
tion is given by the expansion: 

53/2W = E;^- (18) 

n— 1 

The main variables in this approach are the condensate 
density po(j) and the phase space distribution function 
/(r, p) of thermal component. The thermal density pT(r) 
can be obtained from /(r,p) by integrating over mo- 
menta. The effective potential T4(r) and the chemical 
potential p are functions of the condensate density and 
of the thermal density. The condensate and thermal den- 
sities can be found iteratively for a given number of atoms 
and condensate fraction by taking into account that the 
total number of atoms is TV = J dr (/3o(r) + pt{y))- 

The SCHFM is known to work well for the homoge- 
neous harmonic trap [36] and for inhomogeneous traps 
with small aspect ratio |37j . In the present work we need 



the SCHFM not only to speed up the preparation of an 
initial state but also because we will use the chemical 
potential and the thermal atoms distribution function in 
the section |IV| to explain some aspects of the studied 
phenomena. 

III. EXPERIMENTAL SYSTEM 

Following [30], we consider here a cloud of Na atoms 
prepared in the 135*1/2,^ = IjTOf = —1) state and 
confined in a harmonic trap with trapping frequencies 
LUx — ijOy — 2t: X 51Hz and oj^ = 27r x 25Hz. The scat- 
tering length for this system is a = 2.75nni. In subse- 
quent calculations, we use the harmonic oscillator length 
£osc — \/hJmLjl = 4.195/xm, time Tqsc = 1/wz = 6.366ms 
and energy Cosc = f^z as space, time and energy units 
(oscillatory units). 

A. Preparation of initial states 

The preparation of an initial state in the harmonic 
trap Vtr(r) = \'m{uj1.x'^ -\- w^y^ -I- (jjIz^) follows the steps 
described in the previous section. An example of such 
state is shown in the first row of Fig|3] In this particular 
case the temperature of the system is lOOnK and the 
condensate fraction is about 20%. We also prepared two 
more initial states corresponding to different condensate 
fractions, temperatures and numbers of atoms, see Table 
|l| When T = 0, the initial state is simply the ground 
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T[nK] 


fcsT[eosc] 


At[eoscl 


250000 


1.0 


0. 


0. 


22.7 


250000 


0.5 


84. 


69.1 


16.2 


250000 


0.2 


100. 


83.7 


12. 



TABLE I. Numerical values of the condensate fraction A^o /N , 
temperature T, thermal energy fcsT, and chemical potential 
p used in our simulations. 

state of the Gross-Pitaevskii equation. 

In the next step, we split the cloud of atoms into two 
approximately equal parts by rising a Gaussian potential 
barrier Vb(r,t) = Vb{v)f{t) at the center of the harmonic 
trap by means of a linear time-ramp /(t), see Fig(2] Such 
a barrier, with height V}, and width 

Vk{r)^Vbe--"'<, (19) 

can be created by optical means using a blue-detuned 
laser light sheet perpendicular to the x-direction. The 
barrier is ramped at a time tg, chosen at the end of the 
initial equilibration phase, and we have fixed the barrier 
rising time at r = 78.54tosc in our simulations. After 
this perturbation, we let the system reach again equi- 
librium in the double-well trap by evolving the state for 
an additional time r. Finally the system is split into 
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FIG. 2. Sketch of the experimental time-sequence. Sohd and 
dashed Hnes correspond to hnear time-ramps f{t) and h{t) 
respectively. 




two separate clouds containing each about 125000 atoms. 
In all our simulations, the barrier parameters are fixed 
at Vb = 432eosc and Wb = 2.529£osc (~ lO.G^im). The 
full width at half-maximum (FWHM) of the barrier is 
Wb = 2Vh[2wb = 4.214sc (« 17.7fj.m). The height of 
the barrier has been chosen much larger than ksT and 
the chemical potential /i so that both thermal and con- 
densed atoms cannot flow through the barrier, see Table 
III 

As we can create the equilibrium state in a double- well 
potential corresponding to different initial temperatures, 
we can also easily prepare our system in a state where 
temperatures in both wells are different. This can be 
done by replacing the zero temperature component in 
the right well by a nonzero temperature cloud as shown 
in the third row of FigjS] The numbers of atoms in each 
well are approximately equal. We have designed all steps 
of the preparation stage of the initial state of two sub- 
systems with different temperatures having in mind a 
possible and realistic experimental realization. Only the 
last step, i.e. replacing the zero temperature component 
in one subsystem by a finite temperature state has to 
be done differently in the experiment. Heating only one 
subsystem localized in a given well could be done by a 
temporal modulation of the well, followed by a thermal- 
ization. 



B. Opening the channel between the two vessels 

Having prepared two subsystems at different temper- 
atures separated by the potential barrier, we can now 
study their dynamics when a thin channel is rapidly 
opened between the two wells. This is done by switch- 
ing on the channel potential Vc{r,t) — Vc{r)h{t), where 
the linear time-ramp h{t) starts after the equilibration of 
the two subsystems created by the barrier, i.e. at time 
to + 2t, see Figj2j Its duration has been fixed to r/lO in 
all our numerical simulations. 

From an experimental point of view, there are various 




FIG. 3. Three-dimensional surface plots of the trapping po- 
tential and averaged atomic column density at the different 
stages of the simulations. The upper row shows the ini- 
tial harmonic trap (left) during the preparation of the initial 
state. The corresponding atomic density at thermal equilib- 
rium is shown on the right side. The second row shows the 
double-well trap obtained by rising the barrier at the center 
of the harmonic trap (left) and the corresponding equilib- 
rium atomic density (right). The third row shows the den- 
sity of atoms in the double-well trap at zero temperature 
(left) and when the temperatures in each well are different 
(right). In this example, the left well contains a pure con- 
densate (T = 0) whereas the condensate fraction in the right 
well is 20% (T = lOOnK). The last row shows the two wells 
connected through a thin channel (left). The corresponding 
atomic density at some stage of the evolution is shown on the 
right. 



ways to create the channel potential Vc(r). For example, 
starting from an harmonic trap, one could use two or- 
thogonal sheets of blue-detuned laser light propagating in 
the {Oy, Oz) plane. These two sheets build together the 
barrier described earlier in this section and by putting 
two obstacles along their direction of propagation, one 
would create two shadows. Their intersection would open 
the desired channel between the two wells but the min- 
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imum channel width would then be constrained by the 
diffraction effects induced by the two obstacles. However 
widths of the order of few /im should be feasible. Alter- 
native methods would be to use TEq.i Hermite-Gaussian 
modes, or properly designed separate traps [STHM] and 
then focus a red-detuned Gaussian beam. The corre- 
sponding channel potential would be: 

,,.2 (y^ + »^) 

where the Rayleigh length xr = fc^w^ of the chan- 
nel laser beam {k^ is the laser wavenumber) has been 
matched to the barrier parameter Wb- For Wc = 3.5/im, 
one would have Wb — 133/im. The sum of the barrier 
potential Vb(r) and of the new channel potential lg(r) 
is shown in the left frame of Fig|4j In this case, the 
opened channel would have two "potholes" separated by 
a relatively small barrier and these spurious wells would 
trap atoms. In order to observe a superfluid flow and the 
fountain effect, one would then have to make sure that 
the chemical potential ^ is larger than this small bar- 
rier height « Vfc/5. We have run numerical simulations 
(not shown here) and checked that the fountain effect is 
indeed present in this case. 

As this spurious trapping would unnecessarily compli- 
cate (but not kill) our proof-of-principle analysis of the 
fountain effect, we have chosen to work with the following 
channel potential in all our numerical simulations: 

T4(r) = -Vb e-(?''+"')/"'c . (21) 

It has the opposite barrier strength Vb, a Gaussian profile 
with width Wc in the (Oy, Oz) plane and same width Wb 
as the barrier potential along Ox. The sum of Vf,(r) and 
Vc{v) creates a smooth channel between the two vessels 
as it is shown in the right frame of Fig|4j 
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X (in units of f„sc) X (in units of 

FIG. 4. Comparison between the combined barrier and chan- 
nel potential obtained by using a focused Gaussian laser beam 
(left frame) and the one used in our simulations (right frame). 
In the first, atoms get trapped in the "potholes" and the num- 
ber of atoms in the vessels has to be increased in order to 
observe the fountain effect. 

The FWHM-width of the channel is = 2^1112^^. 
The final shape of the total potential (harmonic trap in- 



cluded) is shown in the last row of Figjsjon the left, while 
a typical example of the column density of the evolving 
atomic cloud is shown on the right. In our subsequent nu- 
merical simulations we will use different channel widths 
Wc to compare the behavior of the thermal flow to the 
superfluid one. 



IV. NUMERICAL RESULTS 



The main observations of this paper concern the time 
evolution of two dilute atomic clouds at two different tem- 
peratures and initially prepared in two different potential 
wells (vessels). At a certain time, a "trench" is dug in 
the potential barrier separating the two vessels and the 
atoms can flow from one vessel to the other through the 
channel which has been opened. For classical systems 
one would expect a heat transport from the hotter cloud 
to the colder one, followed by a fast thermalization pro- 
cess. The hot vessel is the potential well on the right and 
it contains only 20% of condensed atoms [T = lOOnK). 
The left well is the cold vessel and it initially contains a 
pure condensate (T = 0). In our simulations, we clearly 
see that, shortly after the two vessels are connected, the 
condensate is flowing very fast from the left cold vessel 
to the right hot vessel as shown in FigjS) In the six initial 
frames we clearly see that the atomic density in the right 
hot vessel is increasing significantly while it is decreasing 
in the left cold vessel. During the same time there is no 
visible transfer of thermal atoms from the hot vessel to 
the cold one. Atoms from the cold vessel are rapidly in- 
jected into the hot vessel. This scenario clearly has the 
flavor of the helium fountain experiment where the su- 
perfluid helium is flowing from the colder big vessel to 
the smaller hot vessel through a thin net of capillaries 
and flnally streams through the small hole in the lid to 
form a jet. In our case we do not see a true fountain 
effect but instead some increase of the atomic density 
in the hot vessel. In fact this physical effect could be 
easily observed in an experiment using standard imaging 
techniques. 

One has to note that, in the original helium fountain 
experiment, there is always a very big reservoir of super- 
fluid atoms. Therefore the fountain effect can persist as 
long as the small vessel is heated. In our case the initial 
number of atoms in each wells is the same. The reservoir 
of cold atoms is thus almost emptied very fast. Then the 
situation gets reversed: the right vessel contains more 
cold atoms than the left one and the atomic cloud starts 
to oscillate back and forth between the two vessels. This 
is clearly seen in Fig[5j where frames 6—11 clearly show 
temporal oscillations of the total atomic density between 
the two vessels (left column). 
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FIG. 5. Snapshots of the time evolution of the column 
atomic densities when the cold left vessel (pure condensate, 
T = 0) and the hot right vessel (condensate fraction 20%, 
T = lOOnK) are connected by a channel. The initial num- 
ber of atoms in each vessel is about 125000. The left column 
of the different frames shows the total atomic density, the 
middle column shows the condensate density and the right 
column shows the density of thermal atoms. The channel 
width is Wc = 2.4£osc (lO/im). The time interval between the 
different frames is about 2.5rosc (~ 15.9ms). 



A. Condensate and thermal component 

The above qualitative findings can be quantified. To 
this end we first have to spht the classical field into con- 
densed and thermal components as described in Sec. |IT] 
The evolution of these components is shown in the middle 
and the right panels of Figjs] The flow starts when the 
channel between the two vessels is fully opened, which 
approximately corresponds to the third frame in Figj5] 



Analyzing the condensate part, we see that its initial 
flow is quite turbulent and a series of shock waves ap- 
pears (frames 3-5). This is because the velocity of 
the superfluid component reaches and exceeds the criti- 
cal velocity. As a result thermal atoms are produced in 
the right well (frame 5, left column) and the condensate 
gets fragmented (frame 5, middle column). After this ini- 
tial turbulent evolution, the flow becomes laminar. We 
have checked that these initial effects are signiflcantly re- 
duced when the temperature difference between the two 
subsystems is smaller. 

A quantitative analysis of the dynamics requires an es- 
timation of temperature of both subsystems. In this dy- 
namical nonequilibrium situation, the notion of tempera- 
ture is questionable. However we can use the condensate 
fraction in the left and the right well as an estimate of 
the 'temperature' of both subsystems. To this end, us- 



ing Eq.(lO), we split the relative occupation numbers of 



the one-particle density matrix modes into left and right 
components: 



OO 
OO 



dy gk{x,y,x,y;t) , 
dy gk{x,y,x,y;t) . 



(22) 



This gives, for each vessel, the condensate, the thermal 
cloud and the total relative occupation numbers: Uq (t), 

it) [X = 



W = Ek=i^kit) and nxit) = nff (t) 



L, R). 

We have drawn the above quantities in Figj6] for two 
different initial condensate fractions in the right well, 
20% (T=100nK) for the top frame and 50% for the bot- 
tom frame (T = 84nK). The thin and thick lines corre- 
spond to the condensate and thermal fractions respec- 
tively. The main observations are the following: (i) The 
initial injection of the left condensate at T = into the 
right well lasts about 47rosc (300ms) in the upper frame, 
and about 31rosc (200ms) in the lower frame; (ii) Af- 
ter the initial injection, the condensate fractions in both 
wells oscillate with a small amplitude around a mean 
value - some condensed atoms flow from one well to the 
other; (iii) The thermal components stay almost constant 
in both wells. 

However, a more detailed analysis shows some ini- 
tial increase of the thermal component during the first 
8 — 16tosc (50 — 100ms) in the left well which is followed 
by a very slow flow of the thermal cloud from the hot to 
the cold part of the system. The initial increase of the 
thermal component can be easily explained. First of all, 
the opening of the channel between the two wells is not 
adiabatic and a thermal fraction is excited in the process 
- see the flrst three panels in Figj5] Secondly, the initial 
flow of the condensed component is very fast and turbu- 
lent so it is another source of thermal excitations. Finally 
a small thermal fraction of atoms is initially present in 
the region of the barrier. 

To prove that the thermo-mechanical effect is indeed 
present in our system, we have to show that mechanical 
equilibrium is reached at once whereas thermal equilib- 
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FIG. 6. Time evolution of the condensate and thermal rela- 
tive occupation numbers in the left and right vessels for two 
different initial condensate fractions in the right well. The 
time unit is Tosc = 6.366ms. Top frame: initial right con- 
densate fraction of 20% (T = lOOnK), final channel width 
Wc = 0.96^osc (~ 4/im). Bottom frame: initial right con- 
densate fraction of 50% (T = 83nK), final channel width 
Wc = ^osc (~ 4.2/^m). Condensate relative occupation num- 
bers: no (^) (thin dotted line) and n^(t) (thin dashed line). 
Thermal relative occupation numbers: nfi^{t) (thick dotted 
line) and n!}{t) (thick dashed line). As one can see, after 
some time, the left and right condensate relative occupation 
numbers oscillate around half the total condensate fraction 
no(t)/2 (thin solid line) whereas the thermal fractions stay 
roughly constant. 



FIG. 7. The upper frame shows the time evolution of the con- 
densate (solid line) and of the thermal (dashed line) fractions 
in the left vessel (thin lines) and in the right vessel (thick 
lines). The time unit is Tosc ~ 6.366ms. The initial conden- 
sate fraction in the right vessel is about 50% (T — 83nK) and 
the channel width is Wc = iosc ~ 4.2/im . As one can see, 
the condensate and thermal components in each vessel never 
equilibrate meaning that the system does not reach thermal 
equilibrium. The lower frame shows the local chemical poten- 
tials calculated in the left (thin line) and in the right (thick 
line) vessels. As one can see, the system is able to reach 
rapidly, in about 31rosc (200ms), a state very close to me- 
chanical equilibrium (/ii, ~ Mfl)- The two distinctive features 
of the helium fountain effect are thus recovered: mechanical, 
but not thermal, equilibrium. 



rium is never reached during the considerably long com- 
putation time of our simulations. To this end we first 
compute and compare the relative condensate and ther- 
mal fractions fo'it) and f^{t) in the left {X — L) and 
the right {X = R) vessels: 



/<f(i) 



frit) 



Nx{t) nx{ty 
Nx{t) 



(23) 
(24) 



The upper frame of FigjT] shows these quantities for an 
initial right condensate fraction of 50% (T = 83nK) 
and the thinest channel width considered here, i.e Wc = 
^osc ~ 4.2/im. It is clearly visible that after 157rosc (Is), 



the condensate fraction in the left well is much larger 
than in the right well. This situation will hold obvi- 
ously even longer. Similarly the thermal components in 
both vessels are very different. This signifies that both 
subsystems are not in thermal equilibrium. During the 
evolution, the initial hot cloud in the right vessel always 
remains much hotter then in the left part. 

To show that the system (almost) reaches mechanical 
equilibrium after a short period of time, we have to con- 



sider the chemical potential defined according to ( 16 ) 



^l{r) = gpo(i-) + 2gpT{r) + Vtrir) 



(25) 



At mechanical equilibrium the chemical potential should 
be position-independent. For comparison we choose two 
positions on opposite sides of the barrier located near the 
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maximum of the initial atomic densities in each wells, 
r_R = ix,y,z) and = (— a:,y, z), and we calculate the 
corresponding local chemical potentials fi^ = and 
= fi(rii). There is however, one technical difhculty. 
The condensate and thermal densities are obtained from 
the diagonalization of the column-averaged one-particle 
density matrix. Therefore, in fact we only know the 2D 
densities in the {Ox, Oy) plane for all eigenmodes. To 
estimate the 3D densities, we need to calculate the Oz- 
width of each eigenmode along the channel. For this we 
average the one-particle density matrix Eq.([7]) along Oy 
ending up with column densities in the {Ox, Oz) plane. 
We extract the Oz-width 'w^{x) for each mode along 
the channel as the FWHM of the corresponding column 
densities. The 3D density is then estimated through 
Pk{x,0,0) — p^^ {x , 0) / wl{x) , where p^^ is the column 
density of the k-th mode in the {Ox, Oy) plane. 
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FIG. 8. Time evolution of the relative occupation numbers of 
the condensate tIq '^(t) (thin lines) and of the thermal cloud 
nip^{t) (thick lines) in the left (dotted lines) and in the right 
(dashed lines) vessels. The time unit is Tobc = 6.366ms. The 
initial condensate fraction in the right well is 50% (T = 83nK) 
and the final width of the channel is Wc = i.OEasc (16.8/im). 
As one can see, after a short initial stage, the right and left 
condensate fractions oscillate around a mean value which is 
half the total condensate fraction no{t)/2 (thin solid line). 
The thermal part, after a while, stays roughly constant but, 
as clearly seen, some part of the thermal cloud flows in phase 
with the condensed atoms. 

Having the 3D densities, we can calculate the chemical 
potentials p^ and pn as the average over a few points 
located around x — —4.6 and x = 4.6£osc (19.3/im) re- 
spectively. The time evolution of these chemical poten- 
tials is shown in the lower frame of Fig[7] Although the 
curves look a bit ragged, we nevertheless see that the 
system rapidly reaches a state very close to mechanical 
equilibrium, p^ ~ pr, in about 31tosc (200ms). In fact 
we observe small out-of-phase oscillations of the chemi- 
cal potentials caused by the back-and-forth oscillations 
of the condensed atoms. 

As a main conclusion of the above discussions, we see 
that our system does present all the three distinctive 
features of the helium fountain experiment: (i) the sys- 



tem cannot achieve thermal equilibrium, (ii) it can only 
achieve mechanical equilibrium, and (iii) the component 
which flows through the very narrow channel connecting 
the two vessels at different temperatures does not trans- 
port heat. 



B. Superfluid and normal component 

To show that our system was not reaching thermal 
equilibrium, we had to divide the classical field into a 
condensate and a thermal component. As the condensate 
component corresponds to the dominant eigenvalue of a 
coarse-grained one-particle density matrix, the thermal 
cloud consists of many modes with relatively small occu- 
pation numbers. This coarse-graining procedure splits 
the system into many different modes. On the other 
hand the standard two-fluid model of the helium fountain 
is solely based on the distinction between a superfluid 
and a normal component. For liquid helium, which is a 
strongly interacting system, there is an essential differ- 
ence between the condensate and the superfluid compo- 
nent. This difference is much less pronounced in the case 
of weakly-interacting trapped atomic condensates, but is 
nevertheless noticeable as pointed out in ^38j where the 
macroscopic excitation of a nonzero momentum mode has 
been studied within the classical fields formalism. The 
authors showed that not only the condensate but also 
phonon-like excitations do participate in the frictionless 
flow. Both the condensate part and these phonon modes 
thus contribute to the superfluid. 

As will be shown in this section, this is also the case 
in the system studied here. A careful reader might have 
already noticed that in FigjS] some part of the thermal 
component oscillates together with the condensate. This 
effect is very small for very narrow channels but is becom- 
ing quite pronounced for wider channels. Fig|8]shows the 
dynamics of the relative occupation numbers of the con- 
densate and of the thermal components when the channel 
width is Wc = 4.0£osc (16.8/Ltm), the initial condensate 
fraction in the right vessel being 50%. It is clearly visible 
that a certain amount of excited atoms is flowing in phase 
with the condensate, back and forth from one vessel to 
the other. 

To explain this behavior, we show in Figj9] the time 
evolution of the relative occupation numbers of the flrst 
seven dominant eigenmodes (in the right well) of the one- 
particle density matrix, the largest occupation number 
corresponding to the condensate. Apart from the con- 
densate, the next two modes in the hierarchy exhibit 
very similar, fast and in-phase oscillations and their oc- 
cupation numbers are significantly larger than those of 
the remaining other modes. These two modes, together 
with the condensate, constitute the three largest coherent 
'pieces' of the system. They might be viewed as three in- 
dependent 'condensates' each of them characterized by a 
particular mode structure, occupation number and heal- 
ing length. 
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FIG. 9. Time evolution of the relative occupation numbers of 
the seven dominant modes in the right vessel of the system 
(top frame). The time unit is Tosc. The solid circles, solid 
squares, and solid diamonds correspond respectively to the 
condensate and to the next two highest occupied modes. The 
remaining four thin lines correspond to the next four modes 
of smaller occupation numbers. The bottom frame shows the 
local healing lengths of these modes. The thin horizontal line 
corresponds to the half the width of the channel, Wc/2. The 
initial condensate fraction in the right vessel is about 50% 
(T = 83nK) and the channel width is Wc = 4.04sc (16.8/im). 



To explain why these three modes can flow freely 
from one vessel to the ether, we calculate the local 
healing length for each of these modes, £,k{x,0,0) — 
1 / y^87ra^)fc|a?7o70) where a; is a distance along the chan- 
nel direction, a is the scattering length and pk is the 3D 
density of the mode estimated through the procedure de- 
scribed previously. The healing lengths are shown in the 
bottom frame of Fig|9] The thin horizontal line corre- 
sponds to half the width of the channel, Wc/2. We im- 
mediately see that the modes flowing together with the 
condensate fulfill the condition: 



2 



(26) 



where ^fe is the "typical" healing length of mode k (for 
example, taken at the middle of the channel). As a gen- 
eral rule, we infer that only modes with a typical heal- 
ing length smaller than half the channel width can flow 
freely. These modes, condensate included, form the su- 



perfluid component. Higher modes, having a typical heal- 
ing length larger than Wc/2, cannot fit into the channel 
and cannot flow: they form the normal component. 

The superfluid fraction and superfluid density are re- 
spectively defined as 



Ps{x,y,t) = YXto \i^k{x,y,t)\'' 



(27) 



where ks is the index of the highest occupied one-particle 
density matrix eigenmode fulfilling < Wc/2. Anal- 
ogously one can define corresponding quantities for the 
normal component, i.e. nN{t) and pN{x,y,t). It is more- 
over convenient to split the superfluid and normal frac- 
tions into their left and right components n^'^(t). 
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FIG. 10. Time evolution of the relative occupation num- 
bers of the superfluid (thin lines) and normal (thick lines) 
components in the left (dotted lines) and right (dashed lines) 
vessels. The time unit is Tosc. The system contains initially 
50% {T — 83nK) of condensed atoms in the right poten- 
tial well and the final width of the channel is Wc = 4.0^osc 
(16.8/im). The thin and thick solid lines show half the to- 
tal superfluid and normal fractions respectively. The normal 
fraction flows smoothly and slowly from the hotter vessel to 
the colder one as expected while the superfluid fraction oscil- 
lates back and forth between the two vessels around a mean 
value being half the total superfluid fraction ns{t)/2 (thin 
solid line). The thick solid line represents half the total nor- 
mal fraction (1 — ns(t))/2 which is never reached by the left 
and right normal components during the time scale of the 
simulation. 



These quantities are shown in Fig 10 It can be seen 
that the normal component flows only very slowly from 
the hotter to the colder well as it is expected for the 
superfluid fountain effect. Comparison with FigjS] fur- 
ther shows that the normal component, contrary to 
the thermal one, does not exhibit any temporal oscil- 
lations. In Fig |ll[ we plot the superfluid column density 
ps{x,y,t) (middle column), the normal column density 
PN{x,y,t) (right column), and the total atomic density 
p{x, y, X, y; t) (left column) for the same parameters as in 
Figj5] The left column is identical as in flg. [5] and is put 
here as a reference. One can clearly see that essentially 
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only the superfluid component travels back and forth be- 
tween the two vessels. The normal component remains 
mainly located in the right hotter vessel and its flow to 
the colder left vessel is almost invisible. 

To estimate the rate of flow of the superfluid fraction, 
we wait for the system to reach its oscillatory regime 
and then fit the (damped) oscillations of the superfiuid 
fraction in the left vessel by: 



DIAGONAL SUPERFLUID NON-SUPERFLUID 
PART PART PART 



F{t) = Asin(27ri/t- 



Bt + C\ (28) 



and extract the oscillation frequency v and the oscilla- 
tion amplitude A of the superfluid flow. The maximal 
superfluid flux through the channel is Fs = 2t: AN v. We 
also fit the slow decrease of the normal fraction in the 
left vessel by the linear function G{t) = at + ji. The 
maximal fiux of the normal atoms is then Fj^ — aN. All 
these quantities are collected in Tables In] and 



Wc[osc.u.] 


A 


!/[Hz] 






Fjv[^] 

-'^ L 777,5 J 


0.96 


0.0254 


3.82 


0.0201 


152 


5.03 


1.2 


0.036 


5.27 


0.0277 


298 


6.93 


2.4 


0.06 


12.6 


0.0502 


1188 


12.55 



TABLE II. The relevant coefficients obtained from our fitting 
procedure and the calculated superfluid and normal rates of 
flow. The initial occupation number of the condensate in the 
right well is 20% (T = lOOnK). 





A 


:/[Hz] 


a 


Fs\—\ 

L m-s -1 


Fn\—\ 


1.0 


0.036 


3.8 


0.013 


215 


3.25 


2.0 


0.078 


9.1 


0.027 


1150 


6.75 


4.0 


0.085 


20.1 




2684 





TABLE III. The relevant coefficients obtained from our fitting 
procedure and the calculated superfluid and normal rates of 
flow. The initial occupation number of the condensate in the 
right well is 50% (T = 83nK). 



Note, that the last row of Table |III| does not contain 
any value for the a coefficient nor for the corresponding 
normal flux Fjq. This is because, for wider channels, 
the rate of flow of the normal component is changing 
significantly in time and fitting the decrease by a linear 
function is no longer reasonable. In this case, the fiow is 



fastest at the beginning as it is visible in Fig 10 

We did not include the value of the coefficients _B, C, 
and /3 in the Tables, even if they increase the precision 
of our fitting procedure, as they are essentially irrelevant 
four our considerations. For channel widths Wc < 5£osc, 
the coefficients 7 turns out to be smaller than the statis- 
tical error (7 ~ 0) and are also not included in Tables [ll] 
and |III[ This observation is in agreement with the fact 
that the dynamics takes place in the coUisionless regime 
as mentioned in the Introduction. 

We see that both superfluid and normal flow rates in- 
crease with the channel width. Moreover, the superfluid 




FIG. 11. Snapshots of the time evolution of the total (left), 
superfluid (middle) and normal (right) column densities. The 
initial condensate fractions are 100% (T = 0) in the left ves- 
sel and about 20% (T = lOOnK) in the right vessel. The 
final channel width is Wc = 2.4£osc (10/im). The time inter- 
val between the frames is about 2.5rosc (15.9ms). As clearly 
seen, the superfluid component oscillates back and forth be- 
tween the two vessels while the normal component is essen- 
tially trapped in the hotter right vessel. 



flow rate is in all cases larger by two or three orders of 
magnitude then the normal one. We expect that the nor- 
mal component behaves like a classical fluid. Therefore, 
its flow rate should correspond to the flux of atoms dis- 
tributed initially according to the classical phase space 
distribution as obtained from the SCHFM equations de- 
scribed in the second section. Our SCHFM calculations 
indeed give a value very close to the one obtained from 
the classical fields dynamics. For example the flux of 
thermal atoms for a system initially prepared with 50% 
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of condensed atoms in the hotter vessel and for a final 
channel width Wc ~ 6.0^osc (25.2/xm) is found to be 
f^v ~ 115.4 atoms/ms. The classical field approxima- 
tion gives a similar result « 182.5 atoms/ms. In- 
deed, the very slow transfer of the normal component is 
a phase space distribution effect - a very small fraction of 
thermal atoms have velocities aligned along the channel. 
On the contrary, the superfluid component is built from 
coherent modes. The coherence of these modes extends 
over the entire two vessels and is established on a short 
time scale of about 16tosc (100ms). 



different temperatures exhibits the two main features of 
the superfluid fountain effect: the mechanical cqiiilibrium 
is obtained almost instantly while the thermal equilib- 
rium is never reached. We have further shown that the 
superfluid component of this system is composed of all 
eigenmodes of the one-particle density matrix having a 
sufficiently small healing length that can fit into the com- 
munication channel. The superfluid flow is at least two 
orders of magnitude faster than the flow of the normal 
component. The slow flow of the normal component can 
be understood as a phase space effect. 



V. CONCLUSION 

In conclusion, we have shown that the analog of the 
thermo-mcchanical effect, observed in the celebrated su- 
perfluid helium II fountain, could be also observed with 
present-day experiments using weakly-interacting degen- 
erate trapped alkali gases. We have proposed a realis- 
tic experimental setup based on a standard harmonic 
conflnement potential and analyzed it with the help of 
the classical fields aproximation method. The trapped 
ultracold gas is first split in two subsytems thanks to 
a potential barrier. Each of the two independent sub- 
systems achieve their own thermal equilibrium, the final 
temperature in the two vessels being different. At a later 
time, a communication channel is opened between the 
two vessels, and the atoms are allowed to flow from one 
vessel to the other. We have shown that the transport 
of atoms between the two subsystems prepared at two 
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